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Indirect CP violation in K — > tttt decays plays a central role in constraining the flavor structure of the Standard Model 
(SM) and in the search for new physics. For many years the leading uncertainty in the SM prediction of this phenomenon 
was the one associated with the nonperturbative strong interaction dynamics in this process. Here we present a fully 
controlled lattice QCD calculation of these effects, which are described by the neutral kaon mixing parameter Bk- We use 
a two step HEX smeared clover-improved Wilson action, with four lattice spacings from a«0.054fm to a«0.093fm and 
pion masses at and even below the physical value. Nonperturbative renormalization is performed in the RI-MOM scheme, 
where we find that operator mixing induced by chiral symmetry breaking is very small. Using fully nonperturbative 
continuum running, we obtain our main result 5^^(3.5 GeV)=0.531(6)stat(2)sys- A perturbative 2-loop conversion yields 



MS-NDR 
K 



(2GeV) = 0.564(6)stat(3)sys(6)pT and Bk = 0.773(8)stat(3)sys(8)pT, which is in good agreement with current 



results from fits to experimental data. 
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1. Introduction 

Neutral kaon mixing is responsible for indirect CP- 
violation in K^TTTT decays. This violation is quantified 
by the parameter e, which is related to quark flavor mix- 
ing parameters and the ratio of hadronic matrix elements 



Bk = 



{K°\0 



AS=2l 



K'>) 



(1) 



where 0'^s=2 ^ ^^^^^^ _ 75)d][s7^(l - 75)^] (cf. [J for 
details). The computation of (fTl) has some advantages over 
direct computations of the O matrix element, such 
as the partial cancellation of statistical and systematic 
uncertainties. 

Note that a precise determination of Bk together with 
experimental measurements of e yields important con- 
straints on the unitary triangle parameters (p, 77) . 



2. Lattice Details 

We use the Nf^2+1 and iV/=3, 2 HEX [5] smeared, 
tree-level clover improved Wilson ensembles generated 
for determining light quark masses [3l [4j . Out of the five 



^CPT is research unit UMR 6207 of the CNRS and of the uni- 
versities Aix-Marseille I, Aix-Marseille II and Sud Toulon- Var, and 
is affiliated with the FRUMAM. 



available lattice spacings, we use the four finest covering 
a range of 0.054 fm < a < 0.093 fm (the low momentum 
cutoff at the largest lattice spacing a«0.116fm does not 
allow for a reliable extraction of the mixing coefficients). 
For the A^/=2-|-l ensembles, the pion masses straddle 
the physical value. Our lattices have sizes as large as 
L~6 fm and finite volume corrections to M^^ [5] are below 
0.5% [1]. We also computed the finite volume corrections 
to Bk using the results from We found that these 
effects are even below the 0.3% level and thus fully under 
control. The iV/=3 configurations are used to compute 
the required renormalization constants nonperturbatively 
using the RI-MOM method [71 |H]. 



3. Nonperturbative renormalization 

Due to explicit chiral symmetry breaking of the Wil- 
son action, the parity even part of operator Q^^^^ = 
{V — A){V — A) — Oi mixes with the other dimension 



six operators O2 — VV — AA, O- 



3/4 



SS T PP and 



O5 = TT where V, A^ S, P, T denote vector, axial-vector, 
scalar, pseudoscalar and tensor AS* = 1 bilinears respec- 
tively. We denote their bare matrix elements {K'~'\Oi\K^) 
by Qi. The renormalization pattern is then given by [8] 



Qicn _ ZikQk — Zij{Sjk + Ajk)Qk. 



(2) 



Hence, the renormalization matrix Zij is decomposed into 
Zij, which, analogously to the continuum renormalization, 
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Figure 1: Nonperturbative running obtained from the continuum ex- 
trapolation of |3| assuming 0(03 a) (orange) or 0{a^) (blue) scaling 
(for details see [4j), divided by the same running computed at NLO. 
We observe agreement of the runnings between fi = 1.8 GeV and 
3.5 GeV within a statistical error which grows to 2% at the lower 
end of the range. 



0(p~^) contributions coming from virtual pion exchanges. 
The dominant corrections are then 0((ap)^) discretization 
errors and an 0{p~'^) term attributed to double pion 
exchanges. We use different fit windows as well as fit 
functions, which include either an (ap)^ term or an 
additional term to estimate systematic effects coming 
from this ambiguity. However, all these effects turn out 
to be very small, as both fit functions give compatible 
results. We also remove the small remaining quark mass 
dependence in the same fit. Fig. [2] shows the mixing 
term Af^ prior to and after removing the discretization 
effects at a«0.077fm lattice spacing. The same data at 
aw0.054fm are also shown for comparison. We observe 
that all mixing coefficients are small and tend to zero as 
a is decreased. 



4. Matrix Elements 



only mixes O2/3 ^-nd O4/5 respectively, and a correction 
Ajfc, quantifying the mixing of different Ok due to chiral 
symmetry breaking. Since Oi docs not mix in the con- 
tinuum, the only relevant terms in the above expression 
are Zn and Ai^, for fc = 2, . . . , 5. Because of ([!]), the 
multiplicative renormalization factor for Bk is given by 
Zbk = Zii/Z\. 

We use the nonperturbative running method [H El HO] to 
circumvent the window-problem of RI-MOM [71 [11] and 
allow matching to other schemes such as MS and RGI 
with small perturbative uncertainties. This means that, 
after an essentially flat extrapolation of Zb^ to vanishing 
quark mass for each lattice spacing, we compute the ratio 



i?«^^,,^(M,3.5GeV) 



zW 



J (3. 5 GeV) 



(3) 



on the three finest lattices at different ji between 1.8 GeV 
and 3.5 GeV. This ratio is then extrapolated to the con- 
tinuum limit yielding the nonperturbative running factor 
from ^ to 3.5 GeV, (/i, 3.5 GeV). Fig. [l] shows that 
our continuum extrapolated results for the running agree 
with NLO perturbation theory [121 [13] in the /i-range con- 
sidered. Thus we set Zg^_^(3.5 GeV) = 3.5 GeV) • 
ZWk pil^)^ where /i is chosen such that ^(/i) can be 
safely extracted for all four lattice spacings. As an addi- 
tional improvement, we subtract a contact term from the 
propagators as described in P^IHTB] . 

The mixing coefficients An, are obtained as described in 
[5], where the additional subtraction [TT] 

Asub/ N TOiAifc(a,mi) - m2Aifc(a,m2) 

An- (a,mi,m2) = 4) 

mi — 1712 

is applied. Here, Aik{a,mi) is the mixing coefhcient 
obtained at quark mass m^. This procedure removes 



To obtain the matrix elements relevant for the computa- 
tion of Bk, we use color-random U(l) wall sources at i = 
and t = T/2 [17\ and vary the time slice r of the opera- 
tor insertion between 1 and T — 1. The relevant operator 
insertions are O1...5. The matrix elements Qi are deter- 
mined by performing constant fits of the time-symmetrized 
plateaus in t as shown in Fig. |3] We use three different 
fitting ranges in order to estimate systematic effects due 
to excited states. Combining the results from these fits 
and the Aij. determined before, we can decompose Bj^ 
into the contributions from the individual Qi. The chi- 
ral symmetry breaking contributions of operators 02....,5 
are 0{as) suppressed relative to that of the leading op- 
erator Oi. However, the corresponding matrix elements 
Q2....,b are chirally enhanced and grow relative to Qi as 
the SU(3) chiral limit is approached. Thus it is impor- 
tant to control the subtraction of these chiral symmetry 
breaking contaminations. The good chiral properties of 
our fermion action mitigate this problem since the contri- 
bution of Qi largely dominates. As shown in Fig. [l] it 
is 98.1(1.2)% for aw0.077fm, M^-120MeV and m, very 
close to its physical value. 



5. Extraction of Physical Bk 

We perform a combined chiral and continuum fit to ex- 
tract the renormalized Bk at the physical mass point and 
in the continuum limit. Since we simulate at or below 
the physical light quark masses, we can perform a safe 
interpolation in the quark masses instead of relying on ex- 
trapolation formulas. For this combined fit we choose the 
following functional form 

Bf{i.^ GeV, X, y, a) = Bf{i.b GeV) • f{x, y) + d(a), (5) 
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Figure 2: Pole subtracted mixing term Al^'' as defined in Ml and extrapolated to the chiral limit at asiO.OTTfm lattice spacing (left). 
The next panel (middle) shows the same data with the 0({ap]^ discretization error removed. The same procedure was applied to A™'' at 
aRi0.054fm (right). The dashed vertical bars indicate the corresponding fitting regions and the horizontal line corresponds to the extracted 
Ai4 along with its statistical la error band. Note the long plateaus in which the data agree with the fit. The extracted mixing coefficient 
tends to zero as a is reduced. 



where f{x, y) with x=M^ and y=2M'^ — describes the 
quark mass dependence. The generic form of / is 



/(a;, y) = (1 + aiox + a2oX^ + a^iy 
+aiixy 



a^x , / X 



(6) 
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We use in total five different fit forms: three Taylor fits 
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Figure 3: Plateau for B^'^'^ (see expression |ljl) at aa!0.065fm, with 
JVf.^~245MeV and time extent T/a=6A, symmetrized with respect 
to t=T/2. Different solid lines represent constant fits for different 
fitting ranges and the shaded bands the corresponding statistical 
error determined on 2000 bootstrap samples. All fits to the different 
ranges agree very well within these errors. 



with different powers in x and y, i.e. with a-^=0 and 
fliOiCioi left free as well as either a2o or an free or set 
to zero. Additionally, we apply two SU(2) xPT fits (cf. 
[IHlIin]) with a^=l and aio ~ 020 = as well as /o,/^,aoi 
free and an either left free or kept fixed to zero. All fits 
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Figure 4; Total contribution of individual Qi, AnQi to Bj^ in % for 
aRi0.077fm and M.;r~120MeV with their errors. The contributions 
from operators Q2 to Q5 are small or even compatible with zero. 



have good fit quality and show full agreement with the ex- 
pected SU(2) chiral behavior for the case a^=f . Since the 
ratio ([1]) is tree-level 0{a) improved, d{a) is chosen propor- 
tional to either a^a or a^, as discussed in [4 . We do not 
include terms whose coefficients are compatible with zero 
in our final fits. In addition, we use the expressions given 
in [6^ to correct the data for the remaining small finite 
volume effects. We further apply two different pion mass 
cuts of 380 MeV and 340 MeV (cf. Together with 

two different fit ranges for pion and kaon mass extractions, 
different fit ranges and functions for obtaining the mixing 
terms and three different scales for extracting the renor- 
malized matrix elements, we end up with 5760 different 
analyses. All of our fit results are very precise and com- 
patible with each other. A sample (Taylor) fit with agU 



3 



0.6 
0.58 
0.56 
m 0.54 

IS) 

ji 0.52 
^ 0.5 
0.48 
0.46 



a-0.093 fm 
a-0.077 fm 
a-0.065 fm 
a-0.054 fm 
cont-limit 



















1 


- — 















0.02 0.04 0.06 0.08 0.1 
M/[GeV^] 

(a) 



0.12 0.14 0.16 





56 r 




0.55 - 




0.54 - 








CD 


0.53 - 


LO 




CO 






0.52 - 






CQ 






0.51 - 




0.5 - 




0.49 - 




0.005 0.01 0.015 0.02 0.025 0.03 0.035 

aga[fm] 

Figure 6: Continuum extrapolation of B^{a, 3. 5GeV), as obtained 
from the fit in Fig. [5] 
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Figure 5: Combined continuum extrapolation and interpolation to 
physical quark mass for a typical Taylor fit out of our 5760 fits. 
The filled square represents our result for B^'(3.5GeV), the dashed 
vertical line the corresponding physical mass. As can be seen from 
panel (a), the interpolation in the light quark mass is mild. The 
slope of the interpolation in rus is somewhat steeper (b), but both 
interpolations arc fully under control. 



scale dependence in Fig. [5] shows an essentially flat chiral 
behaviour and extremely small discretization effects. This 
is also evident from the continuum extrapolation shown 
in Fig. [6j Using the method from [U [50] for a controlled 
determination of all systematics as well as the statistical 
error, our full nonpcrturbative main result reads 



B'^\3.5GeY) = 0.5308(56)stat(23) 



sys, 



(7) 



where the individual contributions to the systematic er- 
ror originate from the subtraction of the mixing terms 
(0.0021), excited state uncertainties (0.0007), extrapolat- 
ing to the continuum and interpolating to physical quark 
masses (both 0.0006), as well as the extraction of the renor- 
malization constant (0.0002). Fig. [t] shows the statis- 
tical and systematic error distributions of our values for 
B^^{3.5 GeV). Both distributions are fairly symmetric and 
clearly show that our final result is dominated by statisti- 
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Figure 7: Statistical distribution (left) and the distribution at- 
tributed to systematic uncertainties (right) for B^{3.5GeV). The 
solid vertical line denotes our final value and the outer and inner 
bands our statistical and systematic errors. This figure emphasizes 
the fact that our overall error is completely dominated by the statis- 
tical uncertainties. 



cal uncertainties. 

For the reader's convenience, we convert our main result 
of Q into the MS-NDR scheme and into the RGI value 
Bk- We do so by using the NLO anomalous dimensions 
of [12 [13] and the beta function at the highest available 
loop order f21]. It is notoriously difficult to reliably as- 
sess the truncation error of a perturbative series, particu- 
larly in the 68% probability sense of our systematic error 
treatment. As the NLO contributions to the conversion 
factors are < 2% and NNLO contributions are typically 
much smaller [22], we add a rather conservative 1% trun- 
cation error, which is larger than a variety of perturbative 
estimates that we have tried. Because this truncation er- 
ror does not fall into our fully controlled systematic error 
framework, we list it separately and do not combine it with 



4 



other systematics. We thus obtam 

5MS-NDR(2 GeV) = 0.5644(59),tat(25),y,(56)pT (8) 
Bk = 0.7727(81),tat(34)sy.(77)pT. (9) 
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Figure 8: Comparison of our result (|9| with the value for Bk ob- 
tained by CKMfitter (ICHEP 10 update to vertical line). The 
dark and light bands correspond to CKMfitter's la and 2a confi- 
dence intervals, respectively. The results from different Nf=2 ( \18\ . 
1st) and Ar/=2-|-l ([l9l|24l|25], 2nd to 4th) lattice computations are 
also shown. 

In Fig. [8] we compare our result to Standard Model 
expectations and other recent lattice results. Our result 
is in good agreement with indirect Bk determinations 
from global Standard Model fits of flavor mixing data 
obtained by CKMfitter (ICHEP 10 updates to [23]). It 
is consistent with expectations obtained by UTfit [53] by 
either including all decays (B/f_aii=0. 94(17)) or neglecting 
the semileptonic channels (i3i<-_no-si=0. 88(13)) in the 
fits. Therefore, we find no evidence for new fundamental 
contributions to indirect CP-violation in K-^tttt decays. 
This is in-line with the findings of [27 . Moreover, we 
hope that the high precision of our result will encourage 
our colleagues responsible for the determination of the 
other contributions to epsilon to work on reducing their 
uncertainties. 
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